Analysis of Argo data South of 60S

By Nick Young and Fabio Viera Machado

In [5]:
from netCDF4 import Dataset, num2date
import glob
from mpl_toolkits.basemap import Basemap
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import calendar
from scipy.stats import describe
from scipy import interpolate
from tqdm import tqdm_notebook as tqdm
import oceansdb
In [62]:
def load_netcdf(path = "Argo_South_60"):
    files = glob.glob("data/" + path + "/**/*.nc", recursive=True)

    lats = []
    lons = []
    datetimes = []
    sst = []
    max_depth = []

    for f in tqdm(files):
        d = Dataset(f)
        lat = d.variables["LATITUDE"][:]
        mask = lat < -60
        lon = d.variables["LONGITUDE"][:]
        lats.extend(lat[mask])
        lons.extend(lon[mask])
        juld = d.variables["JULD"][:]
        units = d.variables["JULD"].getncattr('units')
        dates = num2date(juld, units, "standard")
        datetimes.extend(dates[mask])
        try:
            sst.extend(d.variables["TEMP_ADJUSTED"][mask,0])
        except:
            sst.extend(np.full(len(mask), np.nan))
        max_depth.extend(np.nanmax(d.variables["PRES_ADJUSTED"], axis=1))

    lats = np.array(lats)
    lons = np.array(lons)
    datetimes = np.array(datetimes)
    sst = np.array(sst)
    np.save(path + "_lats", lats)
    np.save(path + "_lons", lons)
    np.save(path + "_dts", datetimes)
    np.save(path + "_sst", sst)
    np.save(path + "_depth", max_depth)
    return lats, lons, datetimes, sst, max_depth

def plot(lats, lons, z = [], title = "Argo profiles south of 60S", cbtitle = "Number of points in bin", vmax=None):
    # setup north polar stereographic basemap.
    # The longitude lon_0 is at 6-o'clock, and the
    # latitude circle boundinglat is tangent to the edge
    # of the map at lon_0. Default value of lat_ts
    # (latitude of true scale) is pole.
    fig = plt.figure(figsize=(15,15))
    m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
    m.drawcoastlines()
    m.fillcontinents(color='black',lake_color='aqua')
    # draw parallels and meridians.
    m.drawparallels(np.arange(-60, 0, 20))
    m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
    #m.drawmapboundary(fill_color='aqua')
    plt.title("{} {}".format(len(lats), title))

    x, y = m(lons, lats)
    if len(z) == 0:
        hh, locx, locy = np.histogram2d(x, y, bins=100)
        # Sort the points by density, so that the densest points are plotted last
        z = np.array([hh[np.argmax(a<=locx[1:]),np.argmax(b<=locy[1:])] for a,b in zip(x,y)])
        #print(describe(z))
        idx = z.argsort()
        x, y, z = x[idx], y[idx], z[idx]
    #m.imshow(heatmap, interpolation='bicubic', cmap="jet")
    m.scatter(x, y, c=z, s=5, alpha=1, cmap="jet", vmax=vmax)
    cb = plt.colorbar()
    cb.set_label(cbtitle, rotation=270)
    plt.show()

def plot_time(dts, title, label):
    fig, ax = plt.subplots(1, 1, figsize=(20,10))
    for i, dt in enumerate(dts):
        k = label[i]
        df = pd.DataFrame({k: dt})
        df = df.groupby(by=df[k].dt.date).count()
        df.plot(ax=ax, style='.', markersize=3)
        ax.xaxis.set_major_locator(mdates.YearLocator())
        ax.set_title(title)
        ax.set_xlabel("Time",fontsize=12)
        ax.set_ylabel("Number of profiles",fontsize=12)
        #ax.set_xlabel("")
        ax.legend()
    plt.show()
In [58]:
argo_lats, argo_lons, argo_dts, argo_sst, argo_depth = load_netcdf()
plot(argo_lats, argo_lons, vmax=100)
plot(argo_lats, argo_lons, argo_sst[:len(argo_lats)], title = "Argo SST", cbtitle = "Temperature °C")
plot(argo_lats, argo_lons, argo_depth[:len(argo_lats)], title = "Argo Depth", cbtitle = "Depth (decibar pressure)")

/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:30: UserWarning: Warning: converting a masked element to nan.
/usr/local/lib/python3.6/dist-packages/numpy/core/numeric.py:591: UserWarning: Warning: converting a masked element to nan.
  return array(a, dtype, copy=False, order=order, subok=True)
In [59]:
print(min(argo_dts), max(argo_dts))
plot_time([argo_dts], "Argo float dailystamps", ["Argo"])
2001-12-21 15:53:45 2019-07-04 19:53:14
In [63]:
seal_lats, seal_lons, seal_dts, seal_sst, seal_depth = load_netcdf("seal")
plot(seal_lats, seal_lons, title = "Seal data south of 60S", vmax=100)
plot(seal_lats, seal_lons, seal_sst[:len(seal_lats)], title = "Seal SST", cbtitle = "Temperature °C")
plot(seal_lats, seal_lons, seal_depth[:len(seal_lats)], title = "Seal Depth", cbtitle = "Depth (decibar pressure)")

/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:30: UserWarning: Warning: converting a masked element to nan.
In [64]:
print(min(seal_dts), max(seal_dts))
plot_time([seal_dts], "Seal data dailystamps", ["Elephant seal"])
2004-02-05 22:38:00.000002 2018-01-14 23:00:00.000003
In [41]:
all_lats = np.concatenate((seal_lats, argo_lats))
all_lons = np.concatenate((seal_lons, argo_lons))
plot(all_lats, all_lons, title="Argo + seal data south of 60S", vmax=100)
c:\users\fvie285\appdata\local\programs\python\python37-32\lib\site-packages\ipykernel_launcher.py:45: MatplotlibDeprecationWarning: 
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
In [65]:
plot_time([seal_dts, argo_dts], "Argo + seal data dailystamps", ["Seal", "Argo"])
In [74]:
def load_netcdf_grid(path = "Argo_South_60", resolution = 1, target_var = "temp", with_depth = False, filter_month = None, filter_season = None):
    if with_depth:
        grid = np.full(shape=[15 * resolution, 360 * resolution, 200], fill_value=np.nan)
    else:
        grid = np.full(shape=[15 * resolution, 360 * resolution], fill_value=np.nan)
    files = glob.glob("data/{}/**/*.nc".format(path), recursive=True)

    for f in tqdm(files[:]):
        d = Dataset(f)
        lat = d.variables["LATITUDE"][:]
        mask = (lat < -60) & (lat > -74.5)
        if filter_month:
            juld = d.variables["JULD"][:]
            units = d.variables["JULD"].getncattr('units')
            dates = num2date(juld, units, "standard")
            datemask = np.array([d.month == filter_month for d in dates])
            mask &= datemask
        if filter_season:
            juld = d.variables["JULD"][:]
            units = d.variables["JULD"].getncattr('units')
            dates = num2date(juld, units, "standard")
            if filter_season == "Summer":
                datemask = np.array([d.month >= 10 or d.month <= 3 for d in dates])
            elif filter_season == "Winter":
                datemask = np.array([d.month > 3 and d.month < 10 for d in dates])
            mask &= datemask
        if any(mask):
            lat = np.round((np.abs(lat[mask]) - 60) * resolution).astype(int)
            lon = d.variables["LONGITUDE"][mask]
            lon = np.round((lon + 180) * resolution).astype(int)
            lon[lon == 360 * resolution] = 0
            if with_depth:
                pres = np.round(d.variables["PRES_ADJUSTED"][mask] / 10).astype(int)
                temp = d.variables["TEMP_ADJUSTED"][mask]
                if target_var == "sal":
                    try:
                        sal = d.variables["PSAL_ADJUSTED"][mask]
                    except:
                        print("No salinity for {}".format(f))
                        continue
            else:
                pres = d.variables["PRES_ADJUSTED"][mask]
                temp = d.variables["TEMP_ADJUSTED"][mask, 0]
            for x in np.unique(lon):
                for y in np.unique(lat):
                    if with_depth:
                        ptmask = (lon == x) & (lat == y)
                        for z in np.unique(pres[ptmask]):
                            if z>= 200 or np.ma.is_masked(z):
                                continue
                            depthmask = pres[ptmask] == z
                            if target_var == "temp":
                                mean_for_pt = np.ma.mean(temp[ptmask][depthmask])
                            elif target_var == "sal":
                                mean_for_pt = np.ma.mean(sal[ptmask][depthmask])
                            if not np.isnan(mean_for_pt) and not np.ma.is_masked(mean_for_pt):
                                grid[y, x, z] = np.nanmean((grid[y, x, z], mean_for_pt))
                    else:
                        ptmask = (lon == x) & (lat == y)
                        if target_var == "n_point":
                            v = np.sum(ptmask)
                            if not np.ma.is_masked(v):
                                grid[y, x] = np.nansum((grid[y, x], v))
                        elif target_var == "temp":
                            temps_at_pt = temp[ptmask]
                            v = np.ma.mean(temps_at_pt)
                            if not np.ma.is_masked(v):
                                grid[y, x] = np.nanmean((grid[y, x], v))
                        elif target_var == "depth":
                            pres_at_pt = pres[ptmask]
                            if len(pres_at_pt):
                                v = np.ma.max(pres[ptmask])
                                if not np.ma.is_masked(v):
                                    grid[y, x] = np.nanmax((grid[y, x], v))
    return grid

def plot_grid(grid, resolution = 1, title="Gridded SST mean for argo data", cbtitle="Temperature °C", vmin=None, vmax=None):
    fig = plt.figure(figsize=(15,15))
    m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
    m.drawcoastlines()
    m.fillcontinents(color='black',lake_color='aqua')
    # draw parallels and meridians.
    m.drawparallels(np.arange(-60, 0, 20))
    m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
    #m.drawmapboundary(fill_color='aqua')
    plt.title("{} at {}° resolution".format(title, 1/resolution))

    x = np.arange(-180, 180.1, 1/resolution)
    y = np.arange(-60, -75.1, 1/-resolution)
    x, y = np.meshgrid(x, y)
    x, y = m(x, y)

    m.pcolormesh(x, y, grid, cmap="jet", vmin=vmin, vmax=vmax)
    cb = plt.colorbar()
    cb.set_label(cbtitle, rotation=270)
    plt.show()

def interp_nans(grid, method='linear'):
    x = np.arange(0, grid.shape[1])
    y = np.arange(0, grid.shape[0])
    mask = np.isnan(grid)
    xx, yy = np.meshgrid(x, y)
    x1 = xx[~mask]
    y1 = yy[~mask]
    newarr = grid[~mask]
    interp = interpolate.griddata((x1, y1), newarr.ravel(), (xx, yy), method=method)
    if method == "linear":
        interp = interp_nans(interp, "nearest")
    return interp
In [66]:
resolution = 2
argo_pt_grid = load_netcdf_grid(target_var="n_point", with_depth = False, resolution = resolution)
seal_pt_grid = load_netcdf_grid(path="seal", target_var="n_point", with_depth = False, resolution = resolution)


In [47]:
combined_pt_grid = np.nansum((argo_pt_grid, seal_pt_grid), axis=0)
combined_pt_grid[np.isnan(argo_pt_grid) & np.isnan(seal_pt_grid)] = np.nan
In [48]:
plot_grid(argo_pt_grid, title="Number of points in gridded argo dataset", cbtitle="Number of points", vmax=50, resolution = resolution)
plot_grid(seal_pt_grid, title="Number of points in gridded seal dataset", cbtitle="Number of points", vmax=50, resolution = resolution)
plot_grid(combined_pt_grid, title="Number of points in gridded argo+seal datasets", cbtitle="Number of points", vmax=100, resolution = resolution)
c:\users\fvie285\appdata\local\programs\python\python37-32\lib\site-packages\ipykernel_launcher.py:70: MatplotlibDeprecationWarning: 
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
In [222]:
resolution = 4
In [223]:
argo_temp_grid = load_netcdf_grid(resolution = resolution)
seal_temp_grid = load_netcdf_grid("seal", resolution = resolution)
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:59: RuntimeWarning: Mean of empty slice
In [224]:
combined_temp_grid = np.nanmean((argo_temp_grid, seal_temp_grid), axis=0)
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:1: RuntimeWarning: Mean of empty slice
  """Entry point for launching an IPython kernel.
In [225]:
plot_grid(argo_temp_grid, resolution = resolution)
plot_grid(seal_temp_grid, title = "Gridded SST mean for seal data", resolution = resolution)
plot_grid(combined_temp_grid, title = "Gridded SST mean for argo+seal data", resolution = resolution)
In [226]:
interp = interp_nans(combined_temp_grid)
In [227]:
plot_grid(interp, resolution = resolution, title="Interpolated gridded SST mean for argo + seal data")
In [155]:
argo_max_depth = load_netcdf_grid(target_var = "depth", resolution = resolution)
seal_max_depth = load_netcdf_grid("seal", target_var = "depth", resolution = resolution)
1350
748
In [159]:
plot_grid(argo_max_depth, resolution = resolution, title="Argo gridded max depth", cbtitle="Depth (dbar)")
plot_grid(seal_max_depth, resolution = resolution, title="Seal gridded max depth", cbtitle="Depth (dbar)", vmax=2000)
In [281]:
for month in range(1, 13):
    argo_temp = load_netcdf_grid(filter_month = month, with_depth=True)
    seal_temp = load_netcdf_grid("seal", filter_month = month, with_depth=True)
    combined_temp = np.nanmean((argo_temp, seal_temp), axis=0)
    np.save("argo_temp_grid_{}".format(month), argo_temp)
    np.save("seal_temp_grid_{}".format(month), seal_temp)
    np.save("combined_temp_grid_{}".format(month), combined_temp)
    interp_temp = interp_nans(combined_temp[:,:,0])
    argomean = np.round(np.nanmean(argo_temp), 2)
    sealmean = np.round(np.nanmean(seal_temp), 2)
    combmean = np.round(np.nanmean(combined_temp), 2)
    interpmean = np.round(np.nanmean(interp_temp), 2)
    title = "Argo 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], argomean)
    plot_grid(argo_temp[:,:,0], title=title)
    title = "Seal 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], sealmean)
    plot_grid(seal_temp[:,:,0], title=title)
    title = "Argo+Seal 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], combmean)
    plot_grid(combined_temp[:,:,0], title=title)
    title = "Interpolated Argo+Seal 0-10 decibar mean for {}. Mean across all points = {}".format(calendar.month_abbr[month], interpmean)
    plot_grid(interp_temp, title=title)
/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:4: RuntimeWarning: Mean of empty slice
  after removing the cwd from sys.path.
In [75]:
for season in ["Summer", "Winter"]:
    argo_temp = load_netcdf_grid(filter_season = season, with_depth=True)
    seal_temp = load_netcdf_grid("seal", filter_season = season, with_depth=True)
    combined_temp = np.nanmean((argo_temp, seal_temp), axis=0)
    np.save("argo_temp_grid_{}".format(season), argo_temp)
    np.save("seal_temp_grid_{}".format(season), seal_temp)
    np.save("combined_temp_grid_{}".format(season), combined_temp)
    interp_temp = interp_nans(combined_temp[:,:,0])
    argomean = np.round(np.nanmean(argo_temp), 2)
    sealmean = np.round(np.nanmean(seal_temp), 2)
    combmean = np.round(np.nanmean(combined_temp), 2)
    interpmean = np.round(np.nanmean(interp_temp), 2)
    title = "Argo 0-10 decibar mean for {}. Mean across all points = {}".format(season, argomean)
    plot_grid(argo_temp[:,:,0], title=title)
    title = "Seal 0-10 decibar mean for {}. Mean across all points = {}".format(season, sealmean)
    plot_grid(seal_temp[:,:,0], title=title)
    title = "Argo+Seal 0-10 decibar mean for {}. Mean across all points = {}".format(season, combmean)
    plot_grid(combined_temp[:,:,0], title=title)
    title = "Interpolated Argo+Seal 0-10 decibar mean for {}. Mean across all points = {}".format(season, interpmean)
    plot_grid(interp_temp, title=title)


/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:4: RuntimeWarning: Mean of empty slice
  after removing the cwd from sys.path.


In [8]:
resolution = 2
try:
    argo_temp_grid_withdepth = np.load("argo_temp_grid_withdepth.npy")
    seal_temp_grid_withdepth = np.load("seal_temp_grid_withdepth.npy")
    combined_temp_grid_withdepth = np.load("combined_temp_grid_withdepth.npy")
except:
    argo_temp_grid_withdepth = load_netcdf_grid(with_depth=True, resolution = resolution)
    seal_temp_grid_withdepth = load_netcdf_grid("seal", with_depth=True, resolution = resolution)
    combined_temp_grid_withdepth = np.nanmean((argo_temp_grid_withdepth, seal_temp_grid_withdepth), axis=0)
    np.save("argo_temp_grid_withdepth", argo_temp_grid_withdepth)
    np.save("seal_temp_grid_withdepth", seal_temp_grid_withdepth)
    np.save("combined_temp_grid_withdepth", combined_temp_grid_withdepth)
    
In [91]:
dbars = [30, 50, 100, 150, 199]
vmax = 5
vmin = -3
for db in dbars:
    dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
    plot_grid(argo_temp_grid_withdepth[:,:,db], title="Gridded argo temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
    plot_grid(seal_temp_grid_withdepth[:,:,db], title="Gridded seal temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
    plot_grid(combined_temp_grid_withdepth[:,:,db], title="Gridded seal+argo temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
    interp = interp_nans(combined_temp_grid_withdepth[:,:,db])
    plot_grid(interp, title="Interpolated gridded seal+argo temperature data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
In [7]:
fig = plt.figure(figsize=(15,15))
m = Basemap(projection='spstere',boundinglat=-55,lon_0=180,resolution='l')
#m.drawcoastlines()
#m.fillcontinents(color='black',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-60, 0, 20))
m.drawmeridians(np.arange(-180, 181, 20), labels=[1,1,0,1])
#m.drawmapboundary(fill_color='aqua')
m.etopo()
plt.title("Bathymetry south of -60S")
c:\users\fvie285\appdata\local\programs\python\python37-32\lib\site-packages\ipykernel_launcher.py:2: MatplotlibDeprecationWarning: 
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
  
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[7]:
Text(0.5, 1.0, 'Bathymetry south of -60S')
In [220]:
def plot_cross_slope(grid, title="Temperature at 179° W", lon=-179, cbtitle="Temperature °C", resolution = 1):
    fig, ax = plt.subplots(1, 1, figsize=(20,10))
    y = np.arange(-60, -75, -1/resolution)
    z = np.arange(0, 2000, 10)
    levels = np.arange(np.nanmin(grid), np.nanmax(grid), .01)
    im = ax.contourf(y, z, grid.T, cmap="jet", levels = levels, extend="both")
    ax.set_xlabel("° Latitude")
    ax.set_ylabel("Depth (dbar)")
    ax.set_title("{} at {} resolution".format(title, 1 / resolution))
    cb = fig.colorbar(im)
    cb.set_label(cbtitle, rotation=270)
    
    db = oceansdb.ETOPO()
    y = np.arange(-60, -75, -.01)
    h = -db['topography'].extract(lat=y, lon=lon)["height"]
    #ax.plot(y, h)
    ax.fill_between(y, 5000, h, color='black')
    ax.set_ylim(2000, 0)
    ax.set_xlim(-60, -74.5)
    
    plt.show()

filtered = combined_temp_grid_withdepth[:,2,:]
interp = interp_nans(filtered)
plot_cross_slope(filtered, resolution = resolution)
plot_cross_slope(interp, title="Interpolated Temperature at 179° W", resolution = resolution)
In [76]:
try:
    argo_sal_grid_withdepth = np.load("argo_sal_grid_withdepth.npy")
    seal_sal_grid_withdepth = np.load("seal_sal_grid_withdepth.npy")
    combined_sal_grid_withdepth = np.load("combined_sal_grid_withdepth.npy")
except:
    argo_sal_grid_withdepth = load_netcdf_grid(with_depth=True, target_var="sal", resolution = resolution)
    seal_sal_grid_withdepth = load_netcdf_grid("seal", with_depth=True, target_var="sal", resolution = resolution)
    combined_sal_grid_withdepth = np.nanmean((argo_sal_grid_withdepth, seal_sal_grid_withdepth), axis=0)
    np.save("argo_sal_grid_withdepth", argo_sal_grid_withdepth)
    np.save("seal_sal_grid_withdepth", seal_sal_grid_withdepth)
    np.save("combined_sal_grid_withdepth", combined_sal_grid_withdepth)
In [45]:
fig, ax = plt.subplots(1, 1, figsize=(20,10))
print(argo_sal_grid_withdepth.shape)
seal_counts = np.sum(~np.isnan(seal_sal_grid_withdepth), axis=2)
argo_counts = np.sum(~np.isnan(argo_sal_grid_withdepth), axis=2)
print(np.argwhere(seal_counts + argo_counts == 319))
depth = np.arange(0, 2000, 10)
ax.set_xlabel("Salinity")
ax.set_ylabel("Depth (dbar)")
ax.set_title("Salinity profile")
ax.set_ylim(2000, 0)
ax.plot(argo_sal_grid_withdepth[11, 466, :], depth)
ax.plot(seal_sal_grid_withdepth[11, 466, :], depth)
ax.legend(["argo", "seal"])
(30, 720, 200)
[[ 11 466]]
Out[45]:
<matplotlib.legend.Legend at 0x7fa4bf8cf9e8>
In [83]:
for season in ["Summer", "Winter"]:
    argo_temp = np.load("argo_temp_grid_{}.npy".format(season))
    seal_temp = np.load("seal_temp_grid_{}.npy".format(season))
    fig, ax = plt.subplots(1, 1, figsize=(20,10))
    depth = np.arange(0, 2000, 10)
    ax.set_xlabel("{} temperature °C".format(season))
    ax.set_ylabel("Depth (dbar)")
    ax.set_title("{} temperature °C profile".format(season))
    ax.set_ylim(2000, 0)
    ax.plot(argo_temp[6, 233, :], depth)
    ax.plot(seal_temp[6, 233, :], depth)
    ax.legend(["argo", "seal"])
In [12]:
dbars = [30, 50, 100, 150, 199]
vmax = 35
vmin = 33.5
for db in dbars:
    dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
    plot_grid(argo_sal_grid_withdepth[:,:,db], title="Gridded argo salinity data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
    plot_grid(seal_sal_grid_withdepth[:,:,db], title="Gridded seal salinity data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
    plot_grid(combined_sal_grid_withdepth[:,:,db], title="Gridded seal+argo salinity data " + dbs, resolution=resolution, vmin=vmin, vmax=vmax)
    interp = interp_nans(combined_sal_grid_withdepth[:,:,db])
    plot_grid(interp, title="Interpolated gridded seal+argo salinity data " + dbs, cbtitle="Salinity", resolution=resolution, vmin=vmin, vmax=vmax)
c:\users\fvie285\appdata\local\programs\python\python37-32\lib\site-packages\ipykernel_launcher.py:70: MatplotlibDeprecationWarning: 
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
In [153]:
def dens_smow(T):
    '''
    Function calculates density of standard mean ocean water (pure water) using EOS 1980.
    INPUT: T = temperature [degree C (ITS-90)]
    OUTPUT: dens = density [kg/m^3]
    '''
    a0 = 999.842594;a1 = 6.793952e-2;a2 = -9.095290e-3;a3 = 1.001685e-4;a4 = -1.120083e-6;a5 = 6.536332e-9
    T68 = T * 1.00024
    dens = (a0 + (a1 + (a2 + (a3 + (a4 + a5*T68)*T68)*T68)*T68)*T68)
    return dens

def dens0(S, T):
    '''
    Function calculates density of sea water at atmospheric pressure
    USAGE:  dens0 = dens0(S,T)
    DESCRIPTION:
        Density of Sea Water at atmospheric pressure using
        UNESCO 1983 (EOS 1980) polynomial.
    INPUT:  (all must have same dimensions)
        S = salinity    [psu      (PSS-78)]
        T = temperature [degree C (ITS-90)]
    OUTPUT:
        dens0 = density  [kg/m^3] of salt water with properties S,T,
        P=0 (0 db gauge pressure)
    '''
    assert S.shape == T.shape
    T68 = T * 1.00024
    # UNESCO 1983 eqn(13) p17.
    b0 = 8.24493e-1;b1 = -4.0899e-3;b2 = 7.6438e-5;b3 = -8.2467e-7;b4 = 5.3875e-9
    c0 = -5.72466e-3;c1 = +1.0227e-4;c2 = -1.6546e-6
    d0 = 4.8314e-4
    dens = (dens_smow(T) + (b0+ (b1+(b2+(b3+b4*T68)*T68)*T68)*T68)*S + (c0+(c1+c2*T68)*T68)*S*np.sqrt(S)+ d0*(S**2) )
    return dens
In [166]:
dens = dens0(combined_sal_grid_withdepth, combined_temp_grid_withdepth)
dbars = [30, 50, 100, 150, 199]
vmax = 1028
vmin = 1026.5

for db in dbars:
    dbs = "from {}-{} dbar".format(db*10, (db+1)*10)
    interp = interp_nans(dens[:,:,db])
    print(np.nanmin(interp), np.nanmax(interp))
    plot_grid(interp, title="Interpolated gridded seal+argo density data " + dbs, cbtitle="Density", resolution=resolution, vmin=vmin, vmax=vmax)
1027.0181758808599 1028.1316975356945
1027.1160229324942 1028.0460193495496
1027.0661376067346 1027.902182726604
1027.1902021401745 1027.8602134606635
1027.0167057697147 1027.8479702251195
In [219]:
plot_cross_slope(dens[:,2,:], title="Density at 179° W", cbtitle="Density", resolution=resolution)
plot_cross_slope(interp_nans(dens[:,2,:]), title="Interpolated density at 179° W", cbtitle="Density", resolution=resolution)
In [ ]: